When setting up a domain and forecast points, it can be useful to quickly visulaize and interact with domain and gage information. This vignette illustrates basic functionality for visualizing the WRF Hydro domain and channel files and potential forecast points.
Load the rwrfhydro package.
library("rwrfhydro")Set the path to the directory of WRF Hydro test cases.
fcPath <- '~/wrfHydroTestCases/Fourmile_Creek_testcase_v2.0'VisualizeDomain()When working with a new model domain, it can be highly informative to explore the physical or geospatial nature of the model through visualization. The WRF-Hydro model has gridded geospatial information for both the LSM and hydrologic components, which are potentially at different spatial resolutions. The hydrologic components may also employ a separate file containing vector information, which we do not discuss here.
This section demonstrates a basic functionality for visualizing the gridded geospatial information in R. These functions suitable for medium sized domains. Proper GIS software is better suited for visualizing large domains and more intensive spatial representations of model components.
The same function VisualizeDomain helps visualize both LSM and hydro components. It plots the information over map data which can be retrieved from a variety of sources, with the help of the ggmap package.
To properly display the spatial information, R currently requires the use of the rgdal and sp packages though these are not required for general installation and use of rwrfhydro. We hope to relax this assumption in the near future. An underlying function GetDomainCoordsProj does the work of obtaining and transforming the projection and/or datum of the spatial data to render it a lat-lon WGS84 map.
We begin by taking a look at the spatial information file for the LSM. First we specify the path to the ‘Fourmile_Creek’ test We pick the “geo” file, which holds the geospatial information for the LSM. The following commands plot the topographic height for each cell over a map.
geoFile <- paste0(fcPath,'/run.ChannelRouting/DOMAIN/geo_em_d01.Fourmile1km.nlcd11.nc')
coordsProj <- GetDomainCoordsProj(geoFile)
ClosureGeo <- VisualizeDomain(geoFile, plotVar='HGT_M', plot=FALSE, plotDf=coordsProj)
mapMargin <- .05*c(-1,1)
closureRtnGeo <-
ClosureGeo(zoom=11, pointsize=9,
grad=topo.colors(15), alpha=.6, maptype='terrain',
subsetRange=range(plotDf$value),
xlim=range(plotDf$long)+mapMargin*1.5,
ylim = range(plotDf$lat)+mapMargin)## Warning: bounding box given to google - spatial extent only approximate.
## converting bounding box to center/zoom specification. (experimental)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.040369,-105.442619&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN
The step of getting coordsProj could be done internally to VisualizeDomain but in some cases (e.g. plotting many variables for a big domain) can speed up the visualization.
For illustration purposes, we set plot=FALSE in the call to VisualizeDomain to skip the default plot output. Typically this first, default plot serves as a basis for tweaking the options to the returned function, named ClosureGeo in this example, which help achieve the desired graphic representation.
Note carefully that VisualizeDomain returns a function, such a returned function is termed a closure: a function with data inside it. (Because our style convention is first-letter capitals for functions, ClosureGeo begins with a capital letter.) The arguments to this function allow most of its main attributes to be tweaked when called. Because the graphics in this function are handled using the ggplot2 package, which has a signifcant learning curve, we design these closures to give the user the 90% solution without th need to learn ggplot2. The return values of the closures give the data and the ggplot objects which can befurther manipulated as necessary.
Critically, zoom and pointsize depend on our domain and even the size of the graphics device. To work more directly with the data and the plot object, see the closureRtnGeo object. The returnComponents argument to the closure function gives even more fine grained output. We’ll illustrate working with the return object below.
A similar set of function calls allows visualization of the hydrologic grids. We look at the gridded channel network file. This time we skip the call to GetDomainCoordsProj.
hydroFile <- paste0(fcPath,'/run.ChannelRouting/DOMAIN/Fulldom_hires_hydrofile.Fourmile100m.nc')
ClosureHydro <- VisualizeDomain(hydroFile, plotVar='CHANNELGRID', plot=FALSE)
closureRtnHydro <-
ClosureHydro(zoom=11, pointsize=1,
## can reference the internal plotDf (or other variables internal to the closure)
location=c(lon=mean(plotDf$long), lat=mean(plotDf$lat)), alpha=.2,
grad=c('white','blue'), maptype='terrain')## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.040415,-105.442563&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN
We highlight use of several of the available arguments to the closure. To see the full set of arguments to the closure, simply invoke the print method on the closure: call the function without paraentheses. This will return the entire code and the arguments are shown at the top with their default values (e.g. argument=default). Reading the coode will reveal two variables which actually come from the VisualizeDomain function: bbox and plotDf. (Because the returned function encloses this data, it is refered to as a closure). Several of the arguments to the closure can reference these variables: location, xlim, ylim, subsetRange.
While the above command shows all the hydrologic grid pixels, many of them are missing values (-9999) because they are not on the gridded channel network. We use the subsetRange keyword to limit the range of values represented in the plot. We again use the xlim and ylim keywords to zoom in. We also limit the color scale to blue.
closureRtnHydro <-
ClosureHydro(zoom=14, pointsize=1,
## can reference the internal plotDf (or other variables internal to the closure)
location=c(lon=mean(plotDf$long), lat=mean(plotDf$lat)),
subsetRange=c(0),
grad='blue', maptype='terrain',
xlim=mean(plotDf$long)+mapMargin/2,
ylim =mean(plotDf$lat)+mapMargin/2/1.5)## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.040415,-105.442563&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN
## Warning: Removed 258 rows containing missing values (geom_point).
To illustrate use of the return values of the above closures, we show how to combine the two plots.
names(closureRtnGeo)## [1] "plotDf" "ggObj"
names(closureRtnHydro)## [1] "plotDf" "ggObj"
closureRtnGeo$ggObj +
geom_point(data=closureRtnHydro$plotDf,
aes(x=long, y=lat), size=.5, shape=21, fill='white', color='darkblue') +
ggtitle('Fourmile Creek, CO - Elevation and Stream Channel') The returnComponents argument to the closures will return more fine-grained objects from the plot in caseu those are useful.
VisualizeChanNtwk()Let’s use look at simulated flows to also see that the flow at the basin outlet is zero. Setup the path to a “CHRTOUT” data file.
chrtFile <- paste0(fcPath,'/run.FullRouting/201306010000.CHRTOUT_DOMAIN1')This is the basic function which shows the flow on the network
LocLinkFun<- VisualizeChanNtwk(chrtFile)## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.037668,-105.437332&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN
## Warning: Removed 1 rows containing missing values (geom_rect).
As with VisualizeDomain, VisualizeChanNtwk returns a closure. You can look at the closure function arguments with
args(LocLinkFun)## function (location = c(lon = mean(range(linkDf$lon)), lat = mean(range(linkDf$lat))),
## zoom = 11, source = "google", maptype = "terrain", padPlot = 0.1,
## gaugeZoom = NULL, clickSelect = FALSE, linkShape = 5, gaugeShape = 4)
## NULL
A main issue is knowing/finding the index of a given point. The click option to the closure lets you click and get the inded. Clicking at the outlet point to see that it has (q=) 0 flow and its index is 350. (We can’t show the interaction in a static document, so you have to try it yourself.)
LocLinkFun(click=TRUE)
## Please click on plot to select nearest point to your click...
## Selected point (in cyan on plot) data:
## ind lon lat q
## 350 -105.325675964 40.0182723999 0 Using the original function, any set of valid indices can be excluded from the data. After excluding the lowest point and getting a new LocLinkFun, clicking on the lowest point reveals that that index is 1.
LocLinkFun<-VisualizeChanNtwk(chrtFile, exclude=350)
LocLinkFun(click=TRUE)
##Please click on plot to select nearest point to your click...
##Selected point (in cyan on plot) data:
## ind lon lat q
## 1 -105.327774048 40.0193214417 1.30229723454We want to add the real-life gauge points and find the nearest stream channel grid cells. You can supply gage points in the following format to the VisualizeChanNtwk function and the nearest neighbors on the channel network are returned.
gaugePts <-
list(orodell =data.frame(lon=254.67374999999998408,
lat=40.018666670000001773),
loganMill =data.frame(lon=254.63508330000001934,
lat=40.042027779999997961),
sunshine =data.frame(lon=254.65122220000000652,
lat=40.05761110000000258) )
LocLinkFun <- VisualizeChanNtwk(chrtFile, gaugePts=gaugePts, exc=350, plot=FALSE)## ** orodell **********
## system chanInd q lon lat
## gauge NA NA -105.326250000 40.0186666700
## model 1 0.580254256725 -105.327774048 40.0193214417
## ** loganMill **********
## system chanInd q lon lat
## gauge NA NA -105.364916700 40.0420277800
## model 259 0.467729389668 -105.364463806 40.0423851013
## ** sunshine **********
## system chanInd q lon lat
## gauge NA NA -105.348777800 40.0576111000
## model 260 0.00161933468189 -105.350837708 40.0423851013
Increase the accuracy of the lon/lat ouput and make the plot that was suppressed in the previous call.
LocLinkFun <- VisualizeChanNtwk(chrtFile, gaugePts=gaugePts, exc=350, plot=FALSE, gaugeAccuracy=17)## ** orodell **********
## system chanInd q lon lat
## gauge NA NA -105.32625000000002 40.018666670000002
## model 1 0.58025425672531128 -105.32777404785156 40.019321441650391
## ** loganMill **********
## system chanInd q lon lat
## gauge NA NA -105.36491669999998 40.042027779999998
## model 259 0.46772938966751099 -105.36446380615234 40.042385101318359
## ** sunshine **********
## system chanInd q lon lat
## gauge NA NA -105.34877779999999 40.057611100000003
## model 260 0.0016193346818909049 -105.35083770751953 40.042385101318359
LocLinkFun()## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.038193,-105.438381&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN
## Warning: Removed 1 rows containing missing values (geom_rect).
Change the amount of padding around the domain and the shape of the gauge symbols.
LocLinkFun(pad=.3, gaugeShape=16)## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.038193,-105.438381&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN
## Warning: Removed 1 rows containing missing values (geom_rect).
Zoom to the orodell gauge.
LocLinkFun(zoom=14, gaugeShape=16, gaugeZoom='orodell', pad=15)## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.018994,-105.327012&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN
## Warning: Removed 1 rows containing missing values (geom_rect).
## Warning: Removed 327 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
Zoom to the Logan Mill gauge
LocLinkFun(zoom=15, gaugeShape=16, gaugeZoom='loganMill', pad=15)## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.042206,-105.36469&zoom=15&size=640x640&scale=2&maptype=terrain&language=en-EN
## Warning: Removed 1 rows containing missing values (geom_rect).
## Warning: Removed 331 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
You can also click in the zoomed view.
LocLinkFun(zoom=15, gaugeShape=16, gaugeZoom='loganMill', pad=15, click=TRUE)